# %% [markdown]
# # 机器学习可解释性分析的糖尿病疾病预测

# %% [markdown]
# ##### 在当今社会，糖尿病正以惊人的速度蔓延，成为全球性的重大公共卫生挑战。据国际糖尿病联盟（IDF）统计，全球糖尿病患者数量持续攀升，不仅给患者个人及其家庭带来沉重的生理、心理和经济负担，也对整个医疗保健系统造成了巨大压力。早期准确地预测糖尿病的发生，对于及时采取干预措施、延缓疾病进展、降低并发症风险以及合理分配医疗资源至关重要。
# ##### 传统的糖尿病预测方法多基于临床经验和简单的统计模型，这些方法在处理复杂的多维度数据时存在一定的局限性，难以充分挖掘数据中的潜在信息。随着机器学习技术的飞速发展，其强大的数据分析和模式识别能力为糖尿病预测带来了新的机遇。众多研究表明，机器学习模型在糖尿病预测方面展现出了较高的准确性和预测效能。
# ##### 本项目在基于基本的数据预处理，模型训练外。引入shape库，深入剖析机器学习模型的决策过程，揭示模型所依赖的关键特征和特征之间的交互作用，使模型的预测结果变得透明和可理解。另一方面，可解释性分析有助于发现与糖尿病发生发展密切相关的潜在生物标志物和风险因素，为糖尿病的病因研究和临床干预提供新的思路和方向。
# ##### 本项目数据集来源于Kaggle平台，特征包含高血压、高胆固醇、胆固醇检查、身体质量指数等，预测标签为是否患病-'1'正类表患病，'0'负类表未患病

# %%
#忽视警告
import warnings
warnings.filterwarnings('ignore')

# %% [markdown]
# ## 认识数据

# %%
import pandas as pd

# 使用pandas的read_csv函数读取CSV文件
data = pd.read_csv(r'D:\项目一信贷风险预测\diabetes.csv')

#随机打乱数据集，防止模型学习到固定的顺序
random_seed = 77
data = data.sample(frac=1, random_state=random_seed)

# 显示前几行数据，以确认数据已正确加载

data.head()

# %%
#样本数量
data.shape

# %%
#数据情况
data.info()

# %% [markdown]
# #### 序号|### 数据集基本信息
# | 序号 | 列名（英文） | 列名（中文） | 非空数量 | 数据类型 |
# | --- | --- | --- | --- | --- |
# | 0 | Diabetes_binary | 糖尿病二元指标 | 253680 non-null | float64 |
# | 1 | HighBP | 高血压 | 253680 non-null | float64 |
# | 2 | HighChol | 高胆固醇 | 253680 non-null | float64 |
# | 3 | CholCheck | 胆固醇检查 | 253680 non-null | float64 |
# | 4 | BMI | 身体质量指数 | 253680 non-null | float64 |
# | 5 | Smoker | 吸烟者 | 253680 non-null | float64 |
# | 6 | Stroke | 中风 | 253680 non-null | float64 |
# | 7 | HeartDiseaseorAttack | 心脏病或心脏病发作 | 253680 non-null | float64 |
# | 8 | PhysActivity | 体育活动 | 253680 non-null | float64 |
# | 9 | Fruits | 水果摄入 | 253680 non-null | float64 |
# | 10 | Veggies | 蔬菜摄入 | 253680 non-null | float64 |
# | 11 | HvyAlcoholConsump | 大量饮酒 | 253680 non-null | float64 |
# | 12 | AnyHealthcare | 任何医疗保健 | 253680 non-null | float64 |
# | 13 | NoDocbcCost | 因费用问题未看医生 | 253680 non-null | float64 |
# | 14 | GenHlth | 总体健康状况 | 253680 non-null | float64 |
# | 15 | MentHlth | 心理健康状况 | 253680 non-null | float64 |
# | 16 | PhysHlth | 身体健康状况 | 253680 non-null | float64 |
# | 17 | DiffWalk | 行走困难 | 253680 non-null | float64 |
# | 18 | Sex | 性别 | 253680 non-null | float64 |
# | 19 | Age | 年龄 | 253680 non-null | float64 |
# | 20 | Education | 教育程度 | 253680 non-null | float64 |
# | 21 | Income | 收入 | 253680 non-null | float64 |
# 

# %% [markdown]
# ### 绘制箱线图 查看离散特征对标签的影响

# %%
import seaborn as sns
import matplotlib.pyplot as plt

# 假设 data 是你的数据框，discrete_vars 是离散变量列表
discrete_vars = ['HighBP', 'HighChol', 'CholCheck', 'Smoker', 'Stroke', 'HeartDiseaseorAttack', 'PhysActivity', 'Fruits','Veggies' , 'HvyAlcoholConsump' ,'AnyHealthcare' ,'NoDocbcCost' ,'Sex' ,'Education']  # 替换为实际的离散变量列名
fig, axes = plt.subplots(1, len(discrete_vars), figsize=(15, 5))

for i, var in enumerate(discrete_vars):
    sns.boxplot(data=data, x=var, y='Diabetes_binary', ax=axes[i])
    axes[i].set_title(f'Boxplot of Diabetes_binary by {var}')
    axes[i].set_xlabel(var)

plt.tight_layout()
plt.show()

# %% [markdown]
# ##### 以 “HighBP”（高血压）对应的箱线图为例，箱体位置较低，说明在高血压不同类别中，患糖尿病的比例相对较为集中在较低水平。

# %% [markdown]
# #### 绘制特征热力图 查看连续特征相关性

# %%
# 导入相应包
import numpy as np
import matplotlib as plt
import seaborn as sns
import matplotlib.pyplot as plt

corr_df = data[['BMI','Age','MentHlth','PhysHlth','GenHlth']].corr()
mask = np.array(corr_df)
mask[np.tril_indices_from(mask)] = False  # mask = np.zeros_like(corr) # mask[np.triu_indices_from(mask)] = True
plt.figure(figsize=(20, 10))
sns.heatmap(corr_df, mask=mask, vmax=0.8, square=True, annot=True)

# %% [markdown]
# ##### 从热力图可以看出不同特征之间的相关强度

# %% [markdown]
# ##### 特征交互良好，无特征之间有较强相关性

# %% [markdown]
# ## 数据预处理

# %%
from sklearn.preprocessing import StandardScaler, OneHotEncoder, LabelEncoder # 用于特征标准化、独热编码和数值编码
from sklearn.model_selection import train_test_split # 用于将数据集划分为训练集和测试集


# 分离特征数据和标签数据
X = data.drop(['Diabetes_binary'], axis=1)  # 特征数据
y = data['Diabetes_binary']  # 标签数据

y

# %%
print(X.columns)

# %% [markdown]
# #### 特征处理并划分训练集和测试集

# %%

import pandas as pd
from sklearn.preprocessing import StandardScaler


df = pd.DataFrame(data)

# 定义连续特征
continuous_features = ['BMI', 'MentHlth', 'PhysHlth', 'Age']

# 分离连续特征
X_continuous = df[continuous_features]

# 创建 StandardScaler 对象
scaler = StandardScaler()

# 对连续特征进行标准化
X_continuous_scaled = scaler.fit_transform(X_continuous)

# 将标准化后的连续特征转换回 DataFrame
X_continuous_scaled_df = pd.DataFrame(X_continuous_scaled, columns=continuous_features)

# 合并标准化后的连续特征与原始数据
df_scaled = df.copy()
X[continuous_features] = X_continuous_scaled_df

# 使用 train_test_split 函数按照 8:2 的比例划分数据集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42)


# %%
# 打印处理后的数据集信息
X_train.info()

# %%
# 查看标签数据的分布
y.value_counts()

# %% [markdown]
# ##### 可以看到，数据集很不平衡，需要过采样平衡数据集

# %%
# 对训练集进行过采样

from imblearn.over_sampling import SMOTE

# 使用 SMOTE 进行过采样,并赋值给自己
smote = SMOTE(random_state=42)
X_train, y_train = smote.fit_resample(X_train, y_train)

# %%
y.value_counts()

# %%
# 打印处理后的数据集信息
X_train.info()

# %% [markdown]
# # 建模训练

# %% [markdown]
# ## 逻辑回归

# %%
from sklearn.metrics import roc_auc_score, confusion_matrix, accuracy_score, classification_report  # 用于模型评估
import seaborn as sns
import matplotlib.pyplot as plt # 绘图
from sklearn.linear_model import LogisticRegression  # 逻辑回归模型
print("Logistic Regression:")
lr = LogisticRegression()
lr.fit(X_train, y_train)
y_pred_lr = lr.predict(X_test)

print("Confusion Matrix:\n", confusion_matrix(y_test, y_pred_lr))
print("Classification Report:\n", classification_report(y_test, y_pred_lr))
print("AUC: ", roc_auc_score(y_test, y_pred_lr))

# %% [markdown]
# ### 经过实测，不平衡数据集未过采样auc为0.55，过采样后auc上升到了0.72，数据集质量得到了提升
# ### 接下来网格搜索法调参，以提升模型精度

# %%

from sklearn.model_selection import train_test_split, GridSearchCV  # 用于数据集划分和超参数调优
# 定义逻辑回归模型
lr = LogisticRegression()

# 定义参数网格
param_grid = {
    'C': [0.01, 0.1, 1, 10, 100],  # 正则化强度的倒数
    'penalty': ['l1', 'l2'],       # 正则化类型
    'solver': ['liblinear']        # 优化算法
}

# 初始化网格搜索
grid_search = GridSearchCV(
    estimator=lr,  # 模型
    param_grid=param_grid,  # 参数网格
    cv=5,  # 五折交叉验证
    scoring='roc_auc',  # 使用 AUC 作为评估指标
    n_jobs=-1  # 使用所有可用的CPU核心
)

# 在训练集上执行网格搜索
grid_search.fit(X_train, y_train)

# 输出最佳参数
print("最佳参数组合:", grid_search.best_params_)

# 使用最佳参数训练模型
best_lr = grid_search.best_estimator_

# 在测试集上进行预测
y_pred_lr = best_lr.predict(X_test)

# 输出混淆矩阵
print("Confusion Matrix:\n", confusion_matrix(y_test, y_pred_lr))

# 输出分类报告
print("Classification Report:\n", classification_report(y_test, y_pred_lr))

# 输出 AUC
y_pred_proba = best_lr.predict_proba(X_test)[:, 1]  # 获取正类的概率
print("测试集 AUC: ", roc_auc_score(y_test, y_pred_proba))

# %% [markdown]
# ##### 经过调参，模型精度提升了7个百分点，性能得到了优化

# %% [markdown]
# ##### 由混淆矩阵，可以看出
# ##### precision（精确率）：0 类精确率 0.94，意味着模型预测为 0 类的样本中，94% 确实是 0 类；1 类精确率 0.30，说明模型预测为 1 类的样本中，仅 30% 是真正的 1 类，1 类精确率较低。
# ##### recall（召回率）：0 类召回率 0.71，即实际为 0 类的样本中，模型能正确识别出 71%；1 类召回率 0.73，即实际为 1 类的样本中，模型能正确识别出 73%。
# ##### f1-score（F1 值）：是精确率和召回率的调和平均，0 类 F1 值 0.81，1 类 F1 值 0.42，综合体现模型在各类上的性能，1 类的综合表现较差。
# ##### support（样本数量）：0 类有 43515 个样本，1 类有 7221 个样本。
# ##### accuracy（准确率）：整体模型的预测准确率为 0.71，即正确预测的样本占总样本的 71%。
# ##### macro avg（宏平均）：对各类别指标简单平均，精确率、召回率、F1 值分别为 0.62、0.72、0.62 ，平等看待各类别。
# ##### weighted avg（加权平均）：根据样本数量对指标加权平均，精确率、召回率、F1 值分别为 0.85、0.71、0.75，因 0 类样本多，结果受 0 类影响大。

# %%
# 计算混淆矩阵
cm = confusion_matrix(y_test, y_pred_lr)

# 使用 seaborn 绘制彩色混淆矩阵
plt.figure(figsize=(8, 6))
sns.heatmap(cm, annot=True, fmt="d", cmap="Blues", cbar=False, 
            xticklabels=['Predicted 0', 'Predicted 1'], 
            yticklabels=['Actual 0', 'Actual 1'])
plt.title("Confusion Matrix")
plt.xlabel("Predicted Label")
plt.ylabel("True Label")
plt.show()

# %% [markdown]
# ## SVM支持向量机

# %%
from sklearn.svm import SVC  # 支持向量机分类模型
# 定义 SVM 模型（使用默认参数）
svm = SVC(probability=True, random_state=42)  # 启用概率估计

# 训练模型
svm.fit(X_train, y_train)

# 在测试集上进行预测
y_pred_proba = svm.predict_proba(X_test)[:, 1]  # 获取正类的概率

# 输出测试集 AUC
print("测试集 AUC: ", roc_auc_score(y_test, y_pred_proba))

# %%
# 导入所需的库
import pandas as pd  # 用于数据处理和分析
import numpy as np  # 用于数值计算
from sklearn.preprocessing import StandardScaler, OneHotEncoder, LabelEncoder  # 用于数据预处理
from sklearn.model_selection import train_test_split, GridSearchCV  # 用于数据集划分和超参数调优
from sklearn.linear_model import LinearRegression  # 线性回归模型
from sklearn.linear_model import LogisticRegression  # 逻辑回归模型
from sklearn.svm import SVC  # 支持向量机分类模型
from sklearn.naive_bayes import GaussianNB  # 高斯朴素贝叶斯模型
from sklearn.ensemble import RandomForestClassifier  # 随机森林分类模型
import xgboost as xgb  # XGBoost模型
import lightgbm as lgb  # LightGBM模型
from sklearn.metrics import roc_auc_score, confusion_matrix, accuracy_score, classification_report  # 用于模型评估
import warnings
# 忽略所有警告
warnings.filterwarnings("ignore")
import seaborn as sns
import matplotlib.pyplot as plt # 绘图

# %% [markdown]
# ## 朴素贝叶斯

# %%
from sklearn.naive_bayes import GaussianNB
# 朴素贝叶斯
print("Naive Bayes:")
nb = GaussianNB()
nb.fit(X_train, y_train)
y_pred_nb = nb.predict_proba(X_test)[:, 1]
print("AUC: ", roc_auc_score(y_test, y_pred_nb))

# %% [markdown]
# ## 网格搜索法调参

# %%

from sklearn.naive_bayes import GaussianNB# 定义朴素贝叶斯模型
nb = GaussianNB()

# 定义参数网格
param_grid = {
    'var_smoothing': [1e-9, 1e-8, 1e-7, 1e-6, 1e-5]  # 方差平滑参数
}

# 初始化网格搜索
grid_search = GridSearchCV(
    estimator=nb,  # 模型
    param_grid=param_grid,  # 参数网格
    cv=5,  # 五折交叉验证
    scoring='roc_auc',  # 使用 AUC 作为评估指标
    n_jobs=-1  # 使用所有可用的CPU核心
)

# 在训练集上执行网格搜索
grid_search.fit(X_train, y_train)

# 输出最佳参数
print("最佳参数组合:", grid_search.best_params_)

# 使用最佳参数训练模型
best_nb = grid_search.best_estimator_

# 在测试集上进行预测
y_pred_proba = best_nb.predict_proba(X_test)[:, 1]  # 获取正类的概率

# 输出测试集 AUC
print("测试集 AUC: ", roc_auc_score(y_test, y_pred_proba))

# %% [markdown]
# ##### 调参结果并不理想

# %% [markdown]
# ## 阈值二次调参

# %%

from sklearn.metrics import accuracy_score, roc_auc_score
import numpy as np
import matplotlib.pyplot as plt
from sklearn.naive_bayes import GaussianNB

# 定义 贝叶斯 模型，使用已知的最佳参数
nb = GaussianNB(var_smoothing =  1e-05)#使用最佳参数

# 训练模型
nb.fit(X_train, y_train)

# 在测试集上进行预测
y_pred_proba = nb.predict_proba(X_test)[:, 1]

# 输出测试集 AUC
print("测试集 AUC: ", roc_auc_score(y_test, y_pred_proba))

# 初始化列表用于存储不同阈值下的准确率
accuracies = []
thresholds = np.arange(0.1, 1, 0.01)

for threshold in thresholds:
    y_pred = (y_pred_proba >= threshold).astype(int)
    accuracy = accuracy_score(y_test, y_pred)
    accuracies.append(accuracy)

# 找到最大准确率及其对应的阈值
max_accuracy_index = np.argmax(accuracies)
optimal_threshold = thresholds[max_accuracy_index]
max_accuracy = accuracies[max_accuracy_index]

print(f"最佳阈值: {optimal_threshold}")
print(f"最大准确率: {max_accuracy}")

# 绘制阈值 - 准确率曲线
plt.plot(thresholds, accuracies)
plt.xlabel('Threshold')
plt.ylabel('Accuracy')
plt.title('Threshold - Accuracy Curve')
plt.show()

# %% [markdown]
# ### 通过遍历找到最佳阈值，模型精度上升了五个百分点，性能再次提升

# %% [markdown]
# ## XGBoost

# %%
# 定义 XGBoost 模型（使用默认参数）
xgb_model = xgb.XGBClassifier(random_state=42)

# 训练模型
xgb_model.fit(X_train, y_train)

# 在测试集上进行预测
y_pred_proba = xgb_model.predict_proba(X_test)[:, 1]  # 获取正类的概率

# 输出测试集 AUC
print("测试集 AUC: ", roc_auc_score(y_test, y_pred_proba))

# %%
# 定义 XGBoost 模型
xgb_model = xgb.XGBClassifier(random_state=42)

# 定义参数网格
param_grid = {
    'max_depth': [3, 6, 9],  # 树的最大深度
    'learning_rate': [0.01, 0.1, 0.3],  # 学习率
    'n_estimators': [100, 200, 300],  # 树的数量
    'subsample': [0.8, 1.0],  # 样本采样比例
    'colsample_bytree': [0.8, 1.0]  # 特征采样比例
}

# 初始化网格搜索
grid_search = GridSearchCV(
    estimator=xgb_model,  # 模型
    param_grid=param_grid,  # 参数网格
    cv=5,  # 五折交叉验证
    scoring='roc_auc',  # 使用 AUC 作为评估指标
    n_jobs=1  # 使用所有可用的CPU核心
)

# 在训练集上执行网格搜索
grid_search.fit(X_train, y_train)

# 输出最佳参数
print("最佳参数组合:", grid_search.best_params_)

# 使用最佳参数训练模型
best_xgb = grid_search.best_estimator_

# 在测试集上进行预测
y_pred_proba = best_xgb.predict_proba(X_test)[:, 1]  # 获取正类的概率

# 输出测试集 AUC
print("测试集 AUC: ", roc_auc_score(y_test, y_pred_proba))

# %% [markdown]
# ## LightGBM

# %%
import lightgbm as lgb
from sklearn.model_selection import train_test_split
from sklearn.metrics import roc_auc_score


# 定义 LightGBM 模型（使用默认参数）
lgb_model = lgb.LGBMClassifier(random_state=42)

# 训练模型
lgb_model.fit(X_train, y_train)

# 在测试集上进行预测
y_pred_proba = lgb_model.predict_proba(X_test)[:, 1]  # 获取正类的概率

# 输出测试集 AUC
print("测试集 AUC: ", roc_auc_score(y_test, y_pred_proba))
# y_pred_proba

# %% [markdown]
# ## 随机森林

# %%
# 定义随机森林模型（使用默认参数）
rf = RandomForestClassifier(random_state=42)

# 训练模型
rf.fit(X_train, y_train)

# 在测试集上进行预测
y_pred_proba = rf.predict_proba(X_test)[:, 1]  # 获取正类的概率

# 输出测试集 AUC
print(" AUC: ", roc_auc_score(y_test, y_pred_proba))

# %%
#定义随机森林模型
rf = RandomForestClassifier(random_state=42)

# 定义参数网格
param_grid = {
    'n_estimators': [50, 100, 200],  # 树的数量
    'max_depth': [None, 10, 20],     # 每棵树的最大深度
    'min_samples_split': [2, 5, 10], # 分裂内部节点所需的最小样本数
    'min_samples_leaf': [1, 2, 4],   # 叶子节点所需的最小样本数
    'max_features': ['sqrt', 'log2'] # 每棵树分裂时考虑的最大特征数
}

# 初始化网格搜索
grid_search = GridSearchCV(
    estimator=rf,  # 模型
    param_grid=param_grid,  # 参数网格
    cv=3,  # 三折交叉验证
    scoring='roc_auc',  # 使用 AUC 作为评估指标
    n_jobs=-1  # 使用所有可用的CPU核心
)

# 在训练集上执行网格搜索
grid_search.fit(X_train, y_train)

# 输出最佳参数
print("最佳参数组合:", grid_search.best_params_)

# 使用最佳参数训练模型
best_rf = grid_search.best_estimator_

# 在测试集上进行预测
y_pred_proba = best_rf.predict_proba(X_test)[:, 1]  # 获取正类的概率

# 输出测试集 AUC
print("测试集 AUC: ", roc_auc_score(y_test, y_pred_proba))

# %% [markdown]
# ##### 经过调参，模型精度上升了3个百分点，模型性能得到了提升

# %% [markdown]
# ## **模型初步评估结果**
# 
# | 模型名称               | 测试集 auc   |
# |------------------------|---------------|
# | 逻辑回归               | 0.7972785482764779|
# | svm支持向量机         | 0.7991634615384615 |
# | 朴素贝叶斯           | 0.7612146995586541      |
# | 随机森林               | 0.8018091996572548     |
# | XGBoot                | 0.8252390337698228     |
# | LightGBM              | 0.7998262246687106    |
# 
# 
# ---
# 
# **说明**
#  **auc**：模型预测准确率，值越接近1表示模型性能越好。

# %% [markdown]
# # 尝试两种模型融合方法提高精度

# %% [markdown]
# ## Stacking模型融合

# %%
from sklearn.ensemble import StackingClassifier
from sklearn.linear_model import LogisticRegression  # 逻辑回归模型
from sklearn.ensemble import RandomForestClassifier  # 随机森林分类模型
from xgboost import XGBClassifier #xgb分类
from lightgbm import LGBMClassifier #lgb分类
# 定义基模型 均使用已调最佳参数
base_models = [
    ('lr', LogisticRegression(
        C =  0.1,  # 正则化强度的倒数
        penalty = 'l2',       # 正则化类型
        solver = 'liblinear'  )),  # 逻辑回归
    ('rf', RandomForestClassifier(
        n_estimators =  200,  # 树的数量
        max_depth =  10,     # 每棵树的最大深度
        min_samples_split = 10, # 分裂内部节点所需的最小样本数
        min_samples_leaf = 2,   # 叶子节点所需的最小样本数
        max_features = 'sqrt' # 每棵树分裂时考虑的最大特征数
    )),  # 随机森林
    ('xgb', XGBClassifier(
        max_depth = 6,  # 树的最大深度
        learning_rate = 0.1,  # 学习率
        n_estimators = 100,  # 树的数量
        subsample =  1.0,  # 样本采样比例
        colsample_bytree = 1.0  # 特征采样比例
    ))  # 线性回归
]

# 定义元模型（LGB）
meta_model = LGBMClassifier(random_state=42)

# 创建 Stacking 回归器
stacking_model = StackingClassifier(
    estimators=base_models,  # 基模型
    final_estimator=meta_model,  # 元模型
    cv=5,  # 交叉验证折数
    n_jobs=1  # 使用所有 CPU 核心
)

# 训练 Stacking 模型
stacking_model.fit(X_train, y_train)

# 预测
y_pred_proba = stacking_model.predict_proba(X_test)[:,1]

# 计算 auc
test_auc = roc_auc_score(y_test, y_pred_proba)

# 输出结果
print("Stacking 模型测试集 auc: ",roc_auc_score(y_test,y_pred_proba))

# 打印 Stacking 模型的参数
print("Stacking 模型参数：")
print(stacking_model.get_params())

# %% [markdown]
# ### 模型加权融合

# %%
base_models = [
    ('lr', LogisticRegression(
        C =  0.1,  # 正则化强度的倒数
        penalty = 'l2',       # 正则化类型
        solver = 'liblinear' 
    )),  # 随机森林
    ('rf', RandomForestClassifier(
        n_estimators =  200,  # 树的数量
        max_depth =  10,     # 每棵树的最大深度
        min_samples_split = 10, # 分裂内部节点所需的最小样本数
        min_samples_leaf = 2,   # 叶子节点所需的最小样本数
        max_features = 'sqrt'  # 最优参数
        
    )),  # XGBoost
    ('xgb', XGBClassifier(
        max_depth = 6,  # 树的最大深度
        learning_rate = 0.1,  # 学习率
        n_estimators = 100,  # 树的数量
        subsample =  1.0,  # 样本采样比例
        colsample_bytree = 1.0  # 最优参数
        
    ))  # LightGBM
]

# 训练基模型
for name, model in base_models:
    model.fit(X_train, y_train)
    print(f"{name} 模型训练完成")

# 获取每个模型的预测结果
y_test_pred_lr = base_models[0][1].predict(X_test)  # 逻辑
y_test_pred_rf = base_models[1][1].predict(X_test)  # 随机森林
y_test_pred_xgb = base_models[2][1].predict(X_test)  # XGBoost

# 定义网格搜索范围
weights_lr = np.linspace(0, 1, 11)  # 随机森林权重范围 [0, 0.1, ..., 1]
weights_rf = np.linspace(0, 1, 11)  # XGBoost 权重范围 [0, 0.1, ..., 1]
weights_xgb = np.linspace(0, 1, 11)  # LightGBM 权重范围 [0, 0.1, ..., 1]

# 初始化最佳 RMSE 和最佳权重
best_rmse = float('inf')
best_weights = None

# 网格搜索
for w_rf in weights_lr:
    for w_xgb in weights_rf:
        for w_lgbm in weights_xgb:
            if w_rf + w_xgb + w_lgbm == 1:  # 确保权重之和为 1
                # 加权平均
                y_test_pred_weighted = (w_rf * y_test_pred_lr +
                                        w_xgb * y_test_pred_rf +
                                        w_lgbm * y_test_pred_xgb)
                # 计算 roc
                test_roc = roc_auc_score(y_test, y_test_pred_weighted)
                # 更新最佳权重
                if test_roc < best_rmse:
                    best_roc = test_roc
                    best_weights = (w_rf, w_xgb, w_lgbm)

# 输出最佳权重和 roc
print(f"最佳权重：随机森林={best_weights[0]:.2f}, XGBoost={best_weights[1]:.2f}, LightGBM={best_weights[2]:.2f}")
print(f"最佳 roc: {best_roc:.4f}")

# %% [markdown]
# ### 效果并不是很理想

# %% [markdown]
# ## 结果可视化

# %%
from sklearn.metrics import roc_curve, auc
from sklearn.model_selection import cross_val_predict
models = {
    "Logistic Regression": LogisticRegression(C=0.1,penalty='l2',solver='liblinear'),

    "Naive Bayes": GaussianNB(var_smoothing=1e-5),
    "Random Forest": RandomForestClassifier(
        max_depth=10, max_features='sqrt', min_samples_leaf=2,
        min_samples_split=10, n_estimators=200
    ),
    "XGBoost": xgb.XGBClassifier(
        colsample_bytree=1.0, learning_rate=0.1, max_depth=6,
        n_estimators=100, subsample=1.0
    ),
    "LightGBM": lgb.LGBMClassifier()
       
}

# 绘制验证集的 ROC 曲线
plt.figure(figsize=(12, 6))
plt.subplot(1, 2, 1)  # 左边图：验证集

for name, model in models.items():
    # 使用交叉验证获取验证集的预测概率
    y_pred_val = cross_val_predict(model, X_train, y_train, cv=5, method='predict_proba')[:, 1]
    fpr, tpr, _ = roc_curve(y_train, y_pred_val)
    roc_auc = auc(fpr, tpr)
    plt.plot(fpr, tpr, lw=2, label=f'{name} (AUC = {roc_auc:.3f})')

plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Validation Set ROC Curves')
plt.legend(loc="lower right")

# 绘制测试集的 ROC 曲线
plt.subplot(1, 2, 2)  # 右边图：测试集

for name, model in models.items():
    # 训练模型并在测试集上进行预测
    model.fit(X_train, y_train)
    y_pred_test = model.predict_proba(X_test)[:, 1]
    fpr, tpr, _ = roc_curve(y_test, y_pred_test)
    roc_auc = auc(fpr, tpr)
    plt.plot(fpr, tpr, lw=2, label=f'{name} (AUC = {roc_auc:.2f})')

plt.plot([0, 1], [0, 1], color='navy', lw=2, linestyle='--')
plt.xlim([0.0, 1.0])
plt.ylim([0.0, 1.05])
plt.xlabel('False Positive Rate')
plt.ylabel('True Positive Rate')
plt.title('Test Set ROC Curves')
plt.legend(loc="lower right")

plt.tight_layout()
plt.show()

# %% [markdown]
# ##### 部分模型性能变化： XGBoost ，验证集上 AUC 为 0.940，而测试集上仅为 0.80，性能有明显下降，可能存在过拟合的情况，即模型在训练数据上表现很好，但在新数据上泛化能力不足。Random Forest 也有类似情况。
# ##### 稳定性分析：Logistic Regression 和 Naive Bayes 在验证集和测试集上的 AUC 值变化相对较小，说明这两个模型的稳定性较好，受数据波动的影响较小。

# %%
import numpy as np
import pandas as pd
import statsmodels.api as sm
import matplotlib.pyplot as plt



# 添加常数项（截距）
X_train_const = sm.add_constant(X_train)

# 使用 statsmodels 拟合逻辑回归模型
logit_model = sm.Logit(y_train, X_train_const)
result = logit_model.fit()

# 输出模型摘要
print(result.summary())

# 提取系数和 p 值
params = result.params  # 系数
p_values = result.pvalues  # p 值

# 创建一个 DataFrame 来存储系数和 p 值
coef_df = pd.DataFrame({
    'Feature': X_train_const.columns,
    'Coefficient': params,
    'P-value': p_values
})

# 打印系数和 p 值
print(coef_df)

# 绘制显著性图（系数条形图）
plt.figure(figsize=(10, 6))
plt.barh(coef_df['Feature'], coef_df['Coefficient'], color='skyblue')
plt.axvline(x=0, color='red', linestyle='--', linewidth=1)
plt.xlabel('Coefficient Value')
plt.ylabel('Feature')
plt.title('Logistic Regression Coefficients')
plt.show()

# 绘制显著性图（p 值条形图）
plt.figure(figsize=(10, 6))
plt.barh(coef_df['Feature'], -np.log10(coef_df['P-value']), color='lightgreen')
plt.axvline(x=-np.log10(0.05), color='red', linestyle='--', linewidth=1, label='P-value = 0.05')
plt.xlabel('-log10(P-value)')
plt.ylabel('Feature')
plt.title('Feature Significance (-log10(P-value))')
plt.legend()
plt.show()

# %% [markdown]
# |分析内容|Logistic Regression Coefficients 图|Feature Significance (-log10 (P-value)) 图|
# | ---- | ---- | ---- |
# |显著特征示例|CholCheck、HighBP|Income、Sex|
# |不显著特征示例|部分系数接近 0 的特征|Veggies、Fruits|
# |系数 / 显著性与目标变量关系|正系数表示正相关，负系数表示负相关|-log10 (P-value) 越大，特征对目标变量的影响越有统计学意义|

# %% [markdown]
# # 机器学习可解释性分析

# %% [markdown]
# ### 引入shape库以LighGBM模型进行可解释性分析

# %%

import shap
import lightgbm as lgb

#重新训练一个新的 LightGBM 模型
lgb_best = lgb.LGBMClassifier(random_state=42)
lgb_best.fit(X_train, y_train)

# 使用 shap 进行解释
# 创建一个 SHAP 解释器，使用训练好的 LightGBM 模型
explainer = shap.Explainer(lgb_best)

# 计算 SHAP 值，使用训练集
shap_values = explainer.shap_values(X_train)

# %%
# 可视化所有类别的SHAP 汇总图
shap.summary_plot(shap_values, X_train)

# %% [markdown]
# ##### 从条形图的长度来看，“BMI” 对模型输出的平均影响程度最大，无论是 Class 0 还是 Class 1，其综合影响在所有特征中最为显著。
# ##### 紧随其后的是 “HighBP”（高血压），说明高血压这个特征对模型结果也有较大影响。
# ##### “Age”（年龄）、“GenHlth”（总体健康状况） 等特征的影响程度也相对较高，在模型中发挥着重要作用。而像 “Stroke”（中风）和 “NoDocbcCost”（无医疗费用相关） 等特征的影响程度较小。

# %%
# SHAP 汇总图设置为条形图，可以显示特征重要性
shap.summary_plot(shap_values[1], X_train, plot_type="bar")

# %% [markdown]
# ##### 在该模型中，BMI,HighBP（高血压）、Age和 GenHlth（总体健康状况）对模型结果有显著影响 是影响预测结果的核心因素

# %%
# SHAP 汇总图默认设置为点图，可以反应特征的重要性和影响方向
shap.summary_plot(shap_values[0], X_train,plot_type="dot", max_display=15)

# %% [markdown]
# ##### 由此汇总图可以看出，GenHlth（总体健康），HighBP（高血压）值越小，shape值越大，模型预测为'1'的概率越大。
# ##### 而收入高低shape值有正有负，对模型预测结果没有明显正相关和负相关性。

# %% [markdown]
# 

# %%
shap.dependence_plot('BMI', shap_values[1], X_train, interaction_index="HighBP")

# %% [markdown]
# ##### BMI 值在较低区间SHAP值大多为负，这表明较低的 BMI 值倾向于对模型输出产生负向影响，即会使模型的预测结果不倾向'1'。（实际 BMI 值不会为负，是上方数据处理后的标准化结果）
# ##### BMI 值较低时，蓝色点（低高血压值）相对较多，且对应的 SHAP 值多为负；而在 BMI 值较高的区间，红色点（高高血压值）的分布逐渐增多。

# %%
shap.dependence_plot('GenHlth', shap_values[1], X_train, interaction_index="HighBP")

# %% [markdown]
# ##### 在 GenHlth 值较高的区间，红色点（高胆固醇值）相对较多，且部分红色点对应的 SHAP 值较高，在总体健康状况较好的情况下，高胆固醇特征可能会进一步增强对模型输出的正向影响。

# %%
shap.dependence_plot('BMI', shap_values[1], X_train, interaction_index="HighChol")

# %% [markdown]
# ##### 在 BMI 值较高的区间，红色点分布逐渐增多，说明高胆固醇样本在高 BMI 情况下更为常见。并且在高 BMI 且高胆固醇（红色点集中区域）时，部分 SHAP 值较高，表明高胆固醇与高 BMI 同时存在时，可能会进一步增强对模型输出的正向影响。

# %% [markdown]
# ##### 通过可解释性分析，BMI，HighBP（高血压），HighCol（高胆固醇），Genhith（身体健康）特征对患糖尿病有着重要的影响。且BMI较大时，往往伴随着高血压，高胆固醇。以上因素综合影响着糖尿病疾病的诱发。

# %% [markdown]
# ## 该项目数据集来源于Kaggle平台

# %% [markdown]
# 


